Topological Excitonic Superfluids in Three Dimensions 



Youngseok Kim, 1 E. M. Hankiewicz, 2 and Matthew J. Gilbert 1 

1 Department of Electrical and Computer Engineering, University of Illinois, Urbana, II, 6180^ 
2 Institut fur Theoretische Physik und Astrophysik, 
Universitat Wilrzburg, Am Hubland, 97074 Wiirzburg, Germany 

We study the equilibrium and non-equilibrium properties of topological dipolar intersurface ex- 
citon condensates within time-reversal invariant topological insulators in three spatial dimensions 
without a magnetic field. We elucidate that, in order to correctly identify the proper pairing 
symmetry within the condensate order parameter, the full three-dimensional Hamiltonian must be 
considered. As a corollary, we demonstrate that only particles with similar chirality play a signifi- 
cant role in condensate formation. Furthermore, we find that the intersurface exciton condensation 
is not suppressed by the interconnection of surfaces in three-dimensional topological insulators as 
the intersurface polarizability vanishes in the condensed phase. This eliminates the surface current 
flow leaving only intersurface current flow through the bulk. We conclude by illustrating how the 
excitonic superfluidity may be identified through an examination of the terminal currents above and 
below the condensate critical current. 



Dipolar excitonic superfluidity (DES) has appeared 
in a veritable manifold of systems including microcavi- 
ties^El, cold atom systemd^l and semiconductor quan- 
tum wells^E31 Within condensed matter, emergent ma- 
terials offer the possibility of finding new DESs. To 
this end, spatially segregated monolayers of graphene 
have been both theoreticall}U^L^ and experimental!}^ 
explored for signatures of excitonic superfluidity. While 
signs of interlayer correlation are experimentally ob- 
served, additional fermionic degrees of freedom, or fla- 
vors, screen the strength of the interlayer interaction^ 
making the observation of DES in graphene multilayers 
challenging. 

The adve nt of time-reversal invariant topological in- 
sulators (Tl) 20 * 21 l has brought renewed interest in find- 
ing DES in condensed matter systems. In sufficiently 
thin TI films, it has been proposed that spatially segre- 
gated surface electrons and holes may bind into a topo- 
logical dipolar intersurface exciton superfluid (TDIES). 
To this point, existing approaches to TDIES have con- 
sidered strictly two-dimensional Dirac surface states sep- 
arated by an insulating spaceiJ^HII] Yet the existence 
of a TDIES in three-dimensions is not a foregone con- 
clusion based on two-dimensional surface state analysis. 
The most obvious drawback being that in a 3D TI, each 
of the surfaces is interconnected and there exists no obvi- 
ous mechanism to segregate the electron and hole layers, 
as in other proposed systems. 

In this Letter, we theoretically study the equilibrium 
and non-equilibrium properties of TDIES in 3D TI and 
show that a stable TDIES may be formed. We link 
this stability of TDIES in 3D TI to the fact that inter- 
surface polarizability vanishes in the TDIES phase for- 
bidding quasiparticle recombination via single particle 
mechanisms. Further, we find that in order to obtain 
the proper form of the condensate order parameter, the 
full 3D Hamiltonian must be used. We propose that the 
TDIES phase may be observed via examination of the 
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FIG. 1: (a): Schematic of topological insulator thin-film sys- 
tem under consideration. The top and bottom surfaces are 
assumed to contain equal numbers of electron and holes, re- 
spectively, (b, c): A plot of order parameter (b) [/A T)Cr and 
(c) £7A r/;Cr as a function of device length at the middle of 
the width. Two matrix elements of A At, At (solid line) and 
A st, st (broken line) are plotted in (b) and two matrix el- 
ements of A^t.st ( s °lid line) and A#t,At (broken line) is 
plotted in (c). The real and imaginary parts are plotted in 
red and blue colors, respectively, and U = 1.5 eV. 



terminal currents via 4-terminal electrical transport mea- 
surements. 

We begin in Fig. [TJa), where we schematically show 
the system we consider. We apply top and bottom po- 
tentials of opposite polarity to induce electrons on the 
top surface and holes on the bottom surface. We at- 
tach contacts CI - C4 to the top and bottom surfaces on 
the left and right sides of the TI through which current 
may be injected and extracted as seen in Fig. [TJa). The 
Hamiltonian for our system is that of a 3D time-reversal 
invariant tP^ 
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with a = x,y,z. In Eq. ([I]), our basis includes two 
orbital components (A, B corresponding for example to 
Pll and P2^ in Bi2Se3 27 ) and the two spin components 
(t, i)- Within our basis, the annihilation operator is 
defined as c k = (c k ^ c k,Bt c k,Aic k ,Bi)- We define the 
requisite gamma matrices in Eq. ([TJ as T a — <j a (g) r x , 
r° = I (g) r z where r a and a a are Pauli matrices for or- 
bital and spin, respectively. Additionally, in Eq. 0, 
d a (k) = hvpka and M(k) = m — (l/2)6/c 2 where i?_f, 
m, and 6 are materials dependent parameters. In Eq. 
([!]), the topological states occur when m/b > 0. 28 In this 
work, we set hvp = 3 eVA, and b = — 9 eVA 2 with an 
applied surface gate bias of V g = 1.0 V which places us 
in the dense electron-hole regime where we expect the 
pairs to form a BCS-type state. In our model the value 
of m is set to —1.5 eV. This value of m ensures that 
the surface states are localized in ^-direction within one 
lattice constant.^ The full single particle Hamiltonian is 
then Fourier transformed into the real space assuming 
low energy excitations to obtain the single particle lat- 
tice Hamiltonian (see Supplementary A) where the lattice 
constant is set to be ao = 3 A. 

With the non-interacting Hamiltonian defined, we now 
specify the intersurface interactions. As long as the 
chemical potential remains within the bulk gap, the sur- 
face state wavefunctions decay exponentially as a func- 
tion of distance from the surface. As such, we may define 
the interactions purely as 2D intersurface interactions be- 
tween the top and bottom surface through a local den- 
sity approximation, Hi nt = —J2<ij>Uijn(i)n(j), with 
n(i) = ^2 s cl(i)c s (i) being the electron density opera- 
tors at a lattice site i with spin and orbital index of 
s = A t, B t, A l,B |. We assume an attractive in- 
tersurface interaction mediated by Coulomb interactions, 
Uij = USij, as such we simplify the intersurface inter- 
action Hamiltonian as, 

H int = -f/]T^ e t(i) es (i)4(i)M*)- (2) 

i s,s' 

In Eq. ([2| , we define annihilation operator of top surface 
(electron layer) as e(i) and bottom surface (hole layer) 
as h(i) at an in-plane lattice site i. Following standard 
mean field decomposition, we may finally obtain our in- 
tersurface interaction contribution as 

H int - U J2J2 [ A s,s'(i)hl,(i)e s (i) 

i s,s> (3) 

+Al s ,(i)et(i)hAi)-\A s ,Ai)\ 2 ~_ ■ 

We define the exciton order parameter a d 29 * 30 * 

A s ,,,(i) = (et(i)Mi)), (4) 
With the order parameter phase expressed as 

<p s , s , = tan" 1 (a» 8 ,/A* 8 ,) , (5) 



where the AJ S , and A y s s , are real and imaginary parts 
of the order parameter A s ^ s >. 

Using the total Hamiltonian of H = Ho + Hi nt , we 
may study the equilibrium properties of the system. We 
turn our focus to the TDIES order parameter which is 
obtained by diagonalizing the system Hamiltonian H 
with the system temperature, T sys = K (see Sup- 
plementary B). As our Hamiltonian has both orbital 
and spin degrees of freedom, there are four possible 
pairings in the order parameter described in Eq. Q. 
To clarify this point, we define exciton order param- 
eter subset as A T5Cr , A r / 5Cr , A T5Cr /, and A r / 5Cr / where 
r (a) stands for the same orbital (spin) pairing, while 
t' (a f ) stands for different orbital (spin) pairing (e.g. 
A r?cr D {A^t,At A Bt)Bt A AiyAi A BiyBi }). We calcu- 
late the order parameter self-consistently in a structure 
of dimensions 99(f) x 33 (y) x 24(i) A, and recognize 
that only two types of the exciton pairing order param- 
eters are non-zero: pairing between the same spin and 
the same orbital (A Tj<T ) shown in Fig. [ljb) and pair- 
ing between the same spin but different orbitals (A r / 5Cr ) 
shown in Fig. [ljc). A r / 5Cr is of particular importance 
as it has not been described in the previous work in- 
volving the effective single surface modeP^. We see that 
A T)Cr is purely real while A r / )Cr is purely imaginary. To 
understand which of these is correct, we calculate the 
ground state energy of system, which is minimized when 
we choose A T ^ a (see Supplementary C). The argument 
of the intersurface phase relationship is also consistent 
with the result. A r / )Cr has y? SjS / = tt/2 whereas A T}CT 
has ip SjS f = 0. Therefore, the purely imaginary order pa- 
rameter is proper as it corresponds to the intersurface 
phase relationship which maximizes intersurface coher- 
ence. With (p s ^ s ' = no TDIES exists as the surfaces 
are completely decoupled. Closer examination of the or- 
der parameter reveals a dependence on quasiparticle chi- 
rality, in which only electrons and holes with identical 
chirality bind. 

Beyond the pairing symmetry, we must understand the 
size of the interaction induced gap in a TDIES. Fig. [2|a) 
shows the size of the self-consistent interaction induced 
gap as a function of the interaction strength, U. Yet we 
know that the intersurface interaction is influenced by 
the dielectric environment. To understand this effect, we 
consider the bulk dielectric constant of a TI (cti) which is 
in contact with top and bottom surface insulating layers 
having dielectric constants of cqi and e<22, respectively 
(see inset of Fig. [2^b)). In this c ase, t he bare intersurface 
Coulomb interaction is given b ^ 24 * 31 !: 

~ , N 87re 2 

where D(q) = (eGi+e T /)(e T /+eG2) e qdj r(e G1 -e T i)(e T i- 
6G2) e~ qd and d is the intersurface separation and q is the 
wavevector. From this, it is possible to estimate interac- 
tion strength in real space, [/^(r), where r is the planar 
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FIG. 2: (a): A plot of the intersurface interaction induced 
energy gap as a function of interaction constant, U with the 
gate bias V g = 1.0 V. The non-zero E g for U = eV orig- 
inates from the finite-size effects, (b): A plot of intersurface 
interaction strength, [/, as a function of dielectric constant of 
topological insulator, cti- The thickness of the TI is fixed at 
d = 24 A. The inset illustrates a schematic dielectric struc- 
ture. 
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FIG. 3: A transmission between CI, C2, C3, and C4 of the 
system (a) without and (b) with intersurface interaction at 
a bias of Vbias — 5 mV . (a): When there is no interaction 
and the system is not gapped, current flows across the device 
resulting in Ti 2 (or T34) dominating the transport, (b): When 
the system is gapped, however, the quasiparticle undergoes a 
similar process to Andreev reflections, resulting in large T13 
and T24 as current flows from top surface to bottom surface 
with no transmission across either of the top T12 or bottom 
surfaces T34. 



radius, using the Fourier transformation of Utb(o)- We 
are particularly interested in the case of r = 0, as Eq. 
([2| only considers local intersurface interactions. Using 
a TI thickness of d = 24 A, we obtain the resultant inter- 
surface interaction strength U t b(0) = 17 as a function of 
the TI dielectric constant, as illustrated in Fig. [2jt>) . In 
this work, we select an intersurface interaction strength 
of U = 1.5 eV in order to ensure a large enough gap 
(Eg ~ 67 meV), that we may observe distinct charac- 
teristics of the condensate phase at a finite-size system 
while not unrealistically large in magnitude. By taking 
the thin- film limit, qd — » 0, we find that Utb(o) g° es to 
Aire 2 1 q(eci + ^2), and is independent of eti- For this 
reason, although the dielectric constant of real material 
is large and results in a reduced intersurface interaction, 
the thin film limit ensures a considerable intersurface in- 
teraction strength. Even in a thin film limit, however, 
the material should have a low level of bulk doping, since 
the intersurface interaction is effectively screened by a 
doped bulk (see Supplementary D). 

With an understanding of the equilibrium properties, 
we now seek to understand the salient non-equilibrium 
properties through the application of the non-equilibrium 
Green's function formalism^. We use a structure of size 
195(f) x 33(y) x 24(5) A in transport calculations to en- 
sure sufficient lateral separation of the surface contacts 
during current injection. We iterate over the Green's 
function and the intersurface interactions until the A SjS / 
reaches self-consistency. Once the self-consistency is 
achieved, the contact and spatially resolved currents (see 
Supplementary E) are calculated. 

One of the key questions concerning the utimate sta- 
bility of the TDIES in 3D arises from the nature of a 
3D TI. In a 3D TI, each of the surfaces is interconnected 
and the single-particle hopping term may easily compete 
with the many-body intersurface interaction. Therefore, 



it may be more energetically favorable for electrons on 
one surface to annihilate holes on the other surface via 
the adjoining surface rather than forming a TDIES. As 
such, elucidating where the current flows in our system 
is one of the most crucial questions. To drive current 
flow, we choose the drag- counter flow bias configuration 
in which Vi = -Vbias, V 2 = Vbias and V 3 = V A = VP 
This configuration will drive an intersurface current flow 
from the top surface to bottom surface on the left side 
of our system and from the bottom surface to the top 
surface on the right hand side. When the system is in a 
TDIES phase, then we expect this to be the only mecha- 
nism for current flow with the superfluid gap forbidding 
transport across either the electron or hole doped surface. 

In Fig. [3J we plot the resultant transmissions from 
each contact to every other contact at an intersurface 
bias which drives a current well below that of the super- 
fluid critical current without and with intersurface in- 
teractions. Here we choose to deal with transmissions 
to track the quasiparticle motion. Normally, one would 
simply examine the terminal currents, however, with the 
surfaces interconnected, there is no way to determine the 
current path to a particular contact. Furthermore, as we 
are well within the linear response regime, the examina- 
tion of individual transmissions will not substantially dif- 
fer within the energy integral and remain a valid method 
to assess current flow. When there are no intersurface 
interactions, as in Fig. [3^a), the transport properties are 
dominated by transmissions directly across the surfaces 
(T12 and T34), however, direct transmissions from the 
top surface to the bottom surface (T13 and T24) and di- 
agonal intersurface transmissions (T14 and T23) are non- 
negligible. This is understood by noting that, although 
we are driving a current across the top surface, the pres- 
ence of gapless states on each of the interconnected sur- 
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FIG. 4: A schematic of the transport directed (x) current 
flow (a) with and (b) without intersurface interaction, (a): 
In the presence of TDIES, current flows between Cl and C3 
on the left and from C2 to C4 on the right, (b): However, in 
the system without gap, all surfaces are interconnected and 
the current flows mainly from Cl to C2. (c): A plot of spa- 
tially resolved transport directed current at the top (left) and 
bottom (right) surfaces with Vbias — 5 mV . The quasiparticle 
tunneling process within the coherence length is manifested as 
an equal amounts of current with a different sign on opposite 
surfaces. The current is normalized by |J ma a:| = 0.11 fiA. 



faces contributes to the total contact current. 

This is to be contrasted with Fig. [3^b) where the in- 
tersurface interactions are included. In this case, direct 
intersurface transmissions dominate the transport char- 
acteristics while transport both across individual surfaces 
and diagonal intersurface transport are negligible. This 
signals the acquisition of a gap corresponding to the for- 
mation of a TDIES. Additionally, we find no diagonal 
intersurface transport yet the side surface connecting the 
top and bottom layer is not gapped. The lack of sur- 
face current flow between the top and bottom surfaces 
lies in the fact that when the TDIES is formed, each of 
the constituent electrons and holes is paired. This forces 
the intersurface polarizability to drop to zero as there 
are no free charges available on either the top or bot- 
tom surfaces to respond to voltage perturbations^. As 
we are in the dense electron-hole regime, we expect the 
intrasurface polarizability to be zero before the onset of 
TDIES. More mathematically, the static intersurface po- 
larization operator^ n = gJ2kdEk n F(Ek) •> where g = 
2 is the number of fermionic degrees of freedom in our 
system, and Ek is the energy of a quasiparticle in the 
condensed state. Since the condensate acquires a gap, 
II must vanish at zero temperature and this is a critical 
insight into the formation of a stable TDIES without the 
necessity of gapping the side surfaces to force intersurface 
segregation. 

In as much as the individual transmissions provide im- 



portant insight into where the currents flow, directional 
current densities provide additional clarity. The process 
for the non-equilibrium conduction is shown in Fig. [4^ a). 
When an electron is injected into the top layer, via Cl, 
within the energy range of the superfluid gap, the elec- 
tron undergoes coherent transport within a character- 
istic distance away from the injecting contact, the co- 
herence length, L c . Beyond L c , the injected electron is 
retrorefiected to the opposite paired surface with oppo- 
site momentum, leading to a significant intersurface cur- 
rent flow, and exiting the system through C3. To con- 
serve current, the system launches the exciton across the 
system which breaks when the exciton reaches the con- 
tacts on the opposite side of the surfaces^^H This is not 
the case in Fig. |4jb) where we do not have a TDIES in 
the system. The top and bottom surfaces are not gapped 
and this allows transfer of charge across the top surface 
from Cl to C2. However, as none of the other surfaces are 
gapped we also expect to see charge transfer from the bot- 
tom contacts reach the top surface contact (e.g. from C4 
to C2). In the Fig. (4^c), we show the spatially resolved 
current in the transport, or x, direction in top and bot- 
tom surfaces with TDIES. In this scenario, we exactly see 
the physical manifestation of the exciton flow described 
in Fig. [4^ a) as an identical amount of current flows in 
top and bottom surfaces with different direction within 
L c w 3 nm away from the contacts. This fact, allows for 
a simple electrical measurement to detect the presence of 
a TDIES. When the system is biased in the drag coun- 
terflow configuration and the system remains below the 
superfluid critical current, the amount of current trans- 
ferred between the two surfaces at each respective side of 
the system is identical (Ii = —Is and —I2 = I4). How- 
ever, above critical current, in condensate gap closes and 
the terminal currents are no longer equal and opposite on 
the respective sides as single particle processes dominate 
the surface and intersurface transport 29 . 

In conclusion, we have examined TDIES in 3D. We find 
that the proper exciton pairing order parameter is purely 
imaginary which is necessary to properly account for the 
system dynamics. Furthermore, we find that the exciton 
order parameter is p-wave and that electrons will only 
bind with holes of the same chirality from different or- 
bit als. We find that TDIES in time-reversal invariant TI 
thin films prefer to bind into excitons with a dominant 
many-body energy which prevents single particle inter- 
surface transport as the intersurface polarizablity drops 
to zero. This allows for the observation of superfluid be- 
havior out of the quantum Hall regime and without the 
need to artificially segregate the surfaces. Finally, we find 
that the presence of a topological superfluid may be elec- 
trically detected by the presence of equal and opposite 
intersurface contact currents. 
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Appendix A: Lattice Model Hamiltonian for 
Time-Reversal Invariant Topological Insulator 

Within a lattice model description, assuming a low 
energy excitation, the d a (k) = hvpk a and M(k) = 
m — (l/2)bk 2 in Eq. (1) is read as 

d a (p) = (Hv F /a )sm(p a a ), 

M(p) = m - 3b/a 2 (Al) 
+ (cos(p x a ) + cosQ^ao) + cos(p z a )), 

for a given lattice constant clq. As a result, Eq. (1) will 
be described as the following form in a lattice model: 



Ho = £4 

k 



36 



m — o ) r° 

at 



( fovp . 



sm(k a a )T a + cos(/c a a )r° 



Cfc. 

(A2) 



The sine and cosine terms in Eq. ( A2 ) are Fourier trans- 
formed into nearest neighbor terms in the real space 
Hamiltonian and, as a result, we obtain 
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Appendix B: Order Parameter Calculation 

It is possible to calculate order parameter A s>s / in 
Eq. (4) via diagonalization of the system Hamiltonian 
H = Ho + Hi n t. Assuming H is hermitian, the cor- 
responding eigenvalue (D) and eigenvector (V) matrix 



satisfy following relationship, 
H = VDV 1 = VDV f 

e 2 •• 



Vi v 2 



vat 



o 



\0 



/( v x t 
( v 2 t 



) 



(Bl) 



where the column vector v m corresponds to a normalized 
eigenvector of the eigenvalue of e m with a total number of 
eigenvalues N. As a result, the eigenstate j m with corre- 
sponding energy e m is connected to the states at a lattice 
site I (ci) via a mapping rule of V as follow, 



m = 


/7i f \ 

72 f 


/( 

( 


V1 t 

v 2 t 


)\ 

) 


/Clt\ 

c 2 t 


= Vt [ C t 






V( 


VjV f 


)/ 


Wv 





(B2) 

Using the matrix identity of VV^ = V = I, it is easy 
to map eigenstates to real space basis as [c^] = V [7^] . 
whose matrix element and its complex conjugate are 



N 



N 



(B3) 



where Vi m stands for a Z,m component matrix element 
of V. As a result, the order parameter in Eq. (4) can be 
calculated as 

A(i, = {c\ cv ) = (jj^ v lml l^ (j^ vLnS^j ^ 

N 

m,m' 
N 

= J2 V l™ V Lf(£m-V), 



(B4) 

where / is the Fermi-Dirac distribution and \i is the bulk 
chemical potential. The first line of the Eq. ( |B4[ ) is from 
Eq. (B3), the second line is from the orthonormality of 
the eigenstates, (7m 7m') 5mm>(rim) = 5mm>f(£m - AO- 
When a lattice site consists of three components, for 
example, I = (x,y,z) and I = (V, y' , z'\ by setting 
i = (x, y) = (x' ', y') with z = t and z' = 6, we can obtain 
the order parameter of A(i) in Eq. (4) at an equilib- 
rium. The calculated order parameter is fed back to the 
intersurface interaction Hamiltonian of Eq. (3), and, as a 
result, the order parameter is obtained self-consistent ly. 
As we point out in the text, our Hamiltonian has both 
orbital and spin degrees of freedom and there are four 
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FIG. 5: The total energy of the system as a function of in- 
teraction constant U. The ground state energy is calculated 
self-consistently at equilibrium and normalized by \Etot\ at 
£7 = eV. Inset magnifies part of the plot. 



possible pairings in the order parameter described in Eq. 
(4). In order to clarify this point, we define exciton or- 
der parameter subset as A T5Cr , A r / 5Cr , A T5Cr /, and A r / 5Cr / 
where r (a) stands for the same orbital (spin) pairing, 
while t' (V) stands for different orbital (spin) pairing 
(e.g. A r ^ D {AA^,At ^5t,Bt &ai,ai A B i,bi})- 



Appendix C: Ground State Energy Calculation 

In the mean-field approach, the intersurface pairing 
scheme whose ground state energy is the lowest will 
be the energetically favorable pairing term. We per- 
form numerical calculations and immediately find that 
the pairing A r?cr / and A r / ?cr / are zero. In order to de- 
termine which pairing term provides the lowest ground 
state among A Tj(J and A r / )Cr , we perform numerical cal- 
culations and obtain the total energy via, 



E tot = (u) = j2 £ <*- u 12 i A ^'(*)i 2 



(CI) 



i,s,s' 



where e a is an eigenvalue of the total Hamiltonian ob- 
tained from Eq. (JbTJ. The index a runs over the occu- 
pied states and A s>s /(z) is defined in Eq. (4). The ground 
state energy in equilibrium is self-consistently calculated 
and the result is presented in Fig. |5j The lowest energy is 
obtained when we choose A r / 5Cr , which is the intersurface 
pairing between the same spins but different orbitals. 



Appendix D: Estimation of the Intersurface 
Interaction U 

The bare Coulomb intersurface interaction strength 
with a consideration of dielectric environment in inset 
of Fig. 2(b) is described by Eq. (6). In order to general- 
ize the argument, we neglect finite-size effect and assume 
x and y as a periodic. Consequently, the real space ex- 
pression is obtained from Fourier transform analysis in 
continuum limit as 



u tb (r) = J d\ u tb ( q y*> 



, 87re 2 ^ (_i)i 

2W ^m^h 1 ^ ' 



(Dl) 



From the first to second line of the Eq. (Dl ), we evaluate 
radial part integration by using the Bessel function^, 



J (s) = 1/(270/^ e*- 



1 = 



(-1)* 



(z) 2 . As we 



=0 2 2l l\ i\ x 

are interested in the on-site intersurface interaction only, 
the in-plane radius is set to be r = and corresponding 
real space expression of the Coulomb interaction in Eq. 
(2) is 



U = U tb (0) 



2?r Jo 



' , 8ne 2 

dq m €TI - 



(D2) 



We use the material parameters of cqi = £gi = 3.9eo 
(Si02), where eo is a vacuum dielectric constant. Assum- 
ing a linear dispersion relationship at the surfaces, the 
Fermi wavevector is calculated as kp = Ep/hvp-, where 
hvjp = 3 eVA and Ep = V g = 1.0 eV (potential induced 
by a gate bias at each surfaces). By setting epi as a 
variational parameter with fixed thickness of d = 24 A, 
the numerical integration is performed and the result is 
shown in Fig. 2(b). 

In addition, we estimate a screening effect of the bulk 
doping on the intersurface interaction. By simplifying the 
problem as a point charge like particle screened by the 
constant background doping as illustrated in Fig. [6](a), 
we calculate Thomas-Fermi wavevectop^ 



Qtf 



2.95 



(r a /a )d/2) 



(D3) 



where the free electron sphere, r s , and the effective Bohr 
radius, ao, are 



/J_\ 1/3 _ 47re TI h 2 
\ 47m / ' a ° m*e 2 



(D4) 



In Eq. (D4), n is electron density and m* is an ef- 
fective mass. Using the effective mass 37 of Bi2Se3 as 
m* ~ 0.155m e and dielectric constant of epi = 100eo, 
the resultant qpF with various doping level is presented 
in Fig. |6|b). The bulk doping effectively screens the 
intersurface interaction as the doping level increases. 



7 



(a) Point charge 




Uniformly 
doped bulk 



10 10 10 
doping density (cm" 5 ) 



FIG. 6: (a): Point charge potential at the surface is screened 
by constant doping in the bulk, (b): The resultant qrF as a 
function of doping level in bulk, x-axis is plotted in a loga- 
rithm scale. 



Appendix E: Terminal and Spatially Resolved 
Current Calculation 



In case of the channel connected to the contact 1 and 2, 
we calculate the current by the multi-channel Landauer- 
Biittiker formula in the limit of coherent transport: 



I(V 12 ) 



2e 
~h 



T 12 {E)[h(E) - f 2 {E)}. (El) 



where V\ 2 is the potential difference between two contact, 
T12 is the transmission and /i( 2 ) is Fermi-Dirac distribu- 
tion of the contact 1 (2). In the non-equilibrium Green's 
function (NEGF) formalism, it is also possible to calcu- 
late the spatially resolved current from point ri to r 2 is 
evaluated by usim j 32 ' 38 1 



I(y x r 2 ) 



ie f dE 

[H(r 12 )(G n (r 2U E) - G p (r 2U E)) 
-tf (r 21 )(G n (r 12 , E) - G p (r 12 , E))} . 



(E2) 



In Eq. (E2), G n ^ p > is the electron (hole) correlation func- 



tions calculated with NEGF method, and G(ri 2 ) and 
H(ri 2 ) represents the off-diagonal block connecting sites 
ri and r 2 which is only nonzero for nearest neighbors. 



* Micro and Nanotechnology Laboratory, University of Illi- 
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